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1.  Introduction 

Air-sea  interaction  has  long  been  recognized  as  an  essential  fac¬ 
tor  in  the  numerical  prediction  of  many  atmospheric  and  oceanic 
phenomena,  including  tropical  cyclones  (TCs).  The  coupling  of 
state-of-the-art  atmosphere,  ocean,  and  wave  models  allows  for  a 
more  precise  representation  of  the  complicated  physical  processes 
occurring  near  the  air-sea  interface.  Fully  air-ocean-wave  coupled 
modeling  has  only  come  to  fruition  recently.  However,  model  cou¬ 
pling,  whether  it  is  air-wave,  air-ocean,  or  ocean-wave  coupling, 
has  long  been  studied  in  research  and  operational  environments. 
For  example,  numerous  European  Centre  for  Medium-Range  Fore¬ 
casts  (ECMWF)  studies  by  Peter  Janssen  and  collaborators  (e.g., 
Janssen  and  Viterbo,  1996;  Janssen  et  al.,  2000,  2002)  present  a 
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global  diagnosis  of  ocean-wave  forecasting,  including  the  effects 
of  ocean  waves  on  the  atmospheric  climate.  Additionally,  three- 
dimensional,  air-ocean-wave  coupled  model  advancements  in 
the  coastal  zone  have  led  to  a  significant  increase  in  understanding 
sediment  transport,  even  in  TC  conditions  (Warner  et  al.,  2008, 
2010). 

The  advancement  of  fully-coupled,  air-ocean-wave  modeling 
has  traditionally  been  hindered  by  uncertainty  in  the  parameteri¬ 
zation  of  air-sea  interactions  and  the  technical  difficulty  of  cou¬ 
pling  complex  modeling  systems.  However,  recent  improvements 
in  the  parameterizations  and  computational  techniques  used  in 
the  Earth  System  Modeling  Framework  (ESMF)  have  allowed  for 
high-resolution  application  of  fully-coupled,  air-ocean-wave  mod¬ 
els.  This  particular  study  utilizes  an  ESMF-based  coupled  system, 
the  Coupled  Ocean  Atmosphere  Mesoscale  Prediction  System- 
Tropical  Cyclone  (COAMPS-TC),  in  which  all  the  interactions  be¬ 
tween  the  atmosphere,  ocean,  and  wave  models  occur  within  an 
ESMF  driver  that  exchanges  information  between  the  three 
models. 
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A  prime  example  of  an  extreme  case  of  air-sea  interaction  oc¬ 
curs  with  the  passage  of  a  TC.  Heat  and  momentum  fluxes,  wind 
stresses,  and  complex  wave-current  interactions  play  a  pivotal  role 
in  both  the  intensity  and  structure  of  TCs.  There  have  been  many 
theoretical  and  numerical  modeling  studies  of  the  ocean  response 
to  TCs.  TC-induced  effects  on  the  ocean  include  SST  cooling,  strong 
surface  currents  (which  can  exceed  2  m  s_1),  upwelling  and  down- 
welling,  large-amplitude  surface  waves,  and  large  internal  gravity 
waves  (Price,  1981;  Ginis  et  al.,  1989;  Jacob  et  al.,  2000;  Morey 
et  al.,  2006;  Fan  et  al.,  2009a, b).  The  atmospheric  response  of  TCs 
to  the  SST  and  upper-ocean  heat  content  is  well  documented  (Knaff 
et.  al.,  2013  and  references  cited  within);  however,  the  atmo¬ 
spheric  response  to  wave  momentum  and  moisture  flux  exchange 
(e.g.,  sea  spray)  is  a  current  focus  of  study  in  air-sea  coupled  mod¬ 
eling.  These  wind-wave-current  interactions  and  coupled  pro¬ 
cesses  in  TCs  are  not  well  understood  and  efforts  are  currently 
ongoing  to  improve  both  the  understanding  and  numerical  model¬ 
ing  of  them  (Drennan  et  al.,  2003;  Hara  and  Belcher,  2004;  Moon 
et  al„  2004a, b;  Chen  et  al.,  2010;  Donelan  et  al.,  2012). 

One  of  the  primary  limiting  factors  in  studies  involving  wind- 
wave-current  interactions  during  TCs  is  insufficient  in  situ  data. 
Although  open-ocean,  tropical  cyclone  wind  and  pressure  mea¬ 
surements  are  frequently  sampled  using  radiosondes,  in  situ  cur¬ 
rent  and  wave  measurements  in  hurricane  conditions  have  been 
sparse  and  lack  sufficient  spatial  and  temporal  resolution  to  inves¬ 
tigate  the  ocean  response  in  detail  (Shay  et  al.,  1989).  Difficulties  in 
forecasting  both  the  initiation  and  track  of  TCs  have  led  researchers 
to  rely  on  chance  encounters  by  ocean  buoys,  moorings,  and 
instruments  that  have  usually  been  deployed  for  other  purposes. 

The  Acoustic  Doppler  Current  Profiler  (ADCP)  array  utilized  in 
this  study  (Fig.  1)  was  part  of  the  Naval  Research  Laboratory’s 
(NRL’s)  Slope  to  Shelf  Energetics  and  Exchange  Dynamics  (SEED) 
project,  an  intensive  measurement  program  of  the  outer  continen¬ 
tal  shelf  and  upper-slope  waters  of  the  northern  Gulf  of  Mexico 
(GOM)  (Mitchell  et  al.  2005;  Wang  et  al.  2005;  Teague  et  al. 
2006).  In  September  2004,  Hurricane  Ivan,  with  maximum 


sustained  winds  greater  than  50ms~\  passed  directly  over  the 
SEED  array.  All  the  moorings  survived  to  provide  one  of  the  best 
data  sets  of  ocean  current  profile  measurements  ever  obtained  in 
such  extreme  conditions.  Historically,  instruments  that  experience 
a  chance  encounter  with  a  major  TC  usually  fail  before  the  most  in¬ 
tense  conditions  are  measured.  For  example,  the  National  Oceanic 
and  Atmospheric  Administration  (NOAA)  National  Weather  Service 
National  Data  Buoy  Center  (NDBC)  buoy  42040  was  directly  im¬ 
pacted  by  the  inner  core  of  Hurricane  Ivan  and  failed  shortly  there¬ 
after.  This  particular  buoy  measured  record  significant  wave 
heights  in  excess  of  15  m  before  the  buoy  stopped  reporting,  read¬ 
ings  that  have  since  been  verified  as  accurate. 

The  in  situ  data  sets  were  also  complemented  by  observations 
from  reconnaissance  aircraft  that  flew  special  missions  into  Hurri¬ 
cane  Ivan  as  the  storm  traversed  the  Caribbean  Sea  and  southern 
GOM.  These  flights  measured  high-resolution  wave  spectra  and 
sea-surface  topography  from  Scanning  Radar  Altimeter  (SRA) 
observations  through  a  joint  effort  between  the  National  Aeronau¬ 
tics  and  Space  Administration  (NASA)  Goddard  Space  Flight  Center 
and  the  NOAA/Atlantic  Oceanographic  and  Meteorological  Labora¬ 
tory/Hurricane  Research  Division  (HRD).  Additionally,  several  sa¬ 
tellite  altimeter  passes  (e.g.,  the  European  Remote  Sensing 
Satellite,  ERS-2)  also  provided  useful  wind-wave  information  for 
the  period  of  interest.  The  combination  of  all  these  observational 
datasets  allows  for  a  unique  analysis  of  the  hurricane’s  wind- 
wave-current  interaction. 

This  study  utilizes  NRL's  COAMPS-TC  (Doyle  et  al.,  2012)  for 
tropical  cyclone  modeling,  while  focusing  on  the  ocean-wave 
model  interactions  during  the  high-wind  TC  conditions  of  Hurri¬ 
cane  Ivan  as  it  traversed  the  GOM.  COAMPS-TC  employs  the 
state-of-the-art  ESMF  to  exchange  wind,  wave,  and  current  infor¬ 
mation  across  the  air-sea  interface  (Campbell  et  al.,  2010).  New 
Simulating  Waves  Nearshore  (SWAN)  wind-to-wave  energy  input 
and  wave-breaking  energy  dissipation  source  functions  (Donelan 
et  al.,  2006;  Babanin  et  al.,  2010;  Rogers  et  al.,  2012)  based  on 
observational  studies  (Young  et  al.,  2005;  Donelan  et  al.,  2006; 
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Fig.  l.  Location  of  the  ADCP  array  in  the  northern  GOM  during  September  2004  (Teague  et  al.,  2007).  Bathymetry  (color  contours),  the  location  of  NDBC  buoy  42040  (red 
triangle),  and  the  best  track  from  the  observations  of  Hurricane  Ivan  (hatched  line)  are  shown.  Inset:  Infrared  satellite  picture  (red  hues  indicate  higher  and  colder  cloud  tops) 
of  Hurricane  Ivan  as  the  TC  approached  the  northern  GOM  coast  on  15  September  2004. 
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Babanin  et  al.,  2007)  are  tested  in  TC  conditions,  as  well  as  a  new 
formulation  of  the  SWAN  ocean  surface  drag  coefficient  derived 
from  high  wind  observations  (Hwang,  2011).  Additionally,  ocean- 
to-wave  model  coupling,  which  includes  the  input  of  the  Navy 
Coastal  Ocean  Model  (NCOM)  ocean  surface  currents  and  water 
levels  to  SWAN,  is  explored  through  observational  comparisons 
of  both  coupled  and  uncoupled  ocean-wave  model  interaction 
simulations  for  Hurricane  Ivan.  Finally,  wave-to-ocean  model  cou¬ 
pling  is  examined  through  the  Stokes  Drift  Current  (SDC)  generated 
by  the  surface  waves  predicted  by  SWAN.  New  idealized  studies 
(e.g.,  Van  Roekel  et  al.,  2012)  show  that  the  angle  between  the  sur¬ 
face  atmospheric  wind  stress  and  the  SDC  is  important  in  the  gen¬ 
eration  of  Langmuir  turbulence  (LT)  and  mixing  in  the  upper  ocean. 
A  misalignment  angle  greater  than  zero  tends  to  decrease  LT  and  a 
misalignment  angle  greater  than  90°  tends  to  inhibit  LT  and  mix¬ 
ing.  In  extreme  situations  such  as  TCs,  the  misalignment  angle  be¬ 
tween  the  SDC  and  the  wind  stress  can  commonly  exceed  90°, 
depending  on  the  direction  of  the  swell  and  the  translation  direc¬ 
tion  of  the  TC  (Sullivan  et  al.,  2012).  The  relationship  between 
the  SDC  and  the  surface  wind  stress  is  examined  for  Hurricane  Ivan 
near  the  eyewall,  where  surface  wind  direction  changes  occur  rap¬ 
idly  over  short  distances. 

Section  2  is  a  general  discussion  of  air-ocean-wave  coupling 
parameters  in  COAMPS-TC.  Section  3  discusses  the  COAMPS-TC 
model  setup  for  Hurricane  Ivan  and  describes  the  in  situ  data.  Sec¬ 
tion  4  presents  an  evaluation  of  the  near-surface  currents  from  the 
Hurricane  Ivan  simulation  based  on  the  ADCP  measurements.  The 
results  of  utilizing  the  new  SWAN  wind  input  and  drag  coefficient 
parameterizations  and  the  NCOM/SWAN  coupling  are  also  pre¬ 
sented  in  Section  4.  Section  5  provides  a  summary  and  a  discussion 
of  future  coupling  applications  with  COAMPS-TC. 

2.  Air-ocean-wave  coupling  methods  in  COAMPS-TC 

A  fully-coupled,  air-ocean-wave  model  can  provide  frequent, 
two-way  feedback  between  the  air,  ocean,  and  wave  components. 
In  COAMPS-TC,  the  atmospheric  component  provides  surface 
atmospheric  pressure  and  momentum,  heat,  and  moisture  fluxes 
to  the  ocean  component,  while  the  ocean  model  returns  its  pre¬ 
dicted  SST  to  the  atmosphere,  which  influences  the  prediction  of 
the  atmospheric  model,  including  the  surface  air-sea  fluxes,  on 
its  next  time  step  (Chen  et  al.,  2010). 

The  feedback  of  the  wave  model  to  the  atmosphere  consists  of  a 
non-dimensional  roughness  length  (Charnock  number),  which  is 
used  to  compute  the  atmospheric  momentum  drag.  Either  of  two 
different  empirical  formulations  for  the  Charnock  number  (Doyle, 
2002;  Moon  et  al.,  2004a, b),  which  are  a  function  of  wave  age 
and  wind  speed,  can  be  selected.  For  this  study,  we  employed  the 
wind-wave  coupling  from  Moon  et  al.  (2004b).  Although  not  a  fo¬ 
cus  of  this  study,  the  ocean-atmosphere  and  wave-atmosphere 
coupling  are  important  to  the  intensity  of  the  TC. 

The  feedback  of  the  ocean  component  to  the  wave  component 
includes  the  input  of  NCOM  surface  currents  and  water  levels  to 
SWAN,  while  the  atmospheric  component  impacts  the  wave  field 
though  the  input  of  the  TC  winds.  Water  levels  can  modify  the 
water  depth  utilized  within  the  wave  model  physics  calculations, 
though  this  effect  is  only  significant  if  the  water  depth  is  suffi¬ 
ciently  shallow  that  the  waves  feel  the  bottom.  More  importantly, 
the  surface  currents  input  to  the  wave  model  alter  the  effective 
wind  speed  (i.e.,  the  wind  speed  relative  to  a  frame  of  reference 
moving  with  the  currents)  and  the  horizontal  shear  of  the  currents 
produces  changes  in  the  length,  height,  and  direction  of  the  waves 
in  a  manner  similar  to  refraction  and  shoaling  by  interaction  with 
variable  bathymetry.  These  are  conservative  processes;  e.g.,  shoal¬ 
ing  is  associated  with  conservation  of  the  wave  energy  flux.  How¬ 
ever,  this  can,  in  turn,  lead  to  non-conservative  effects.  For 


example,  waves  that  become  more  (or  less)  steep  due  to  interac¬ 
tion  with  currents  will  be  more  (or  less)  likely  to  break,  and  this 
breaking,  also  called  whitecapping,  is  a  highly  nonlinear  and 
non-conservative  process.  The  effect  of  whitecapping  is  discussed 
later  in  Section  4,  which  shows  results  for  Hurricane  Ivan  with 
and  without  ocean-to-wave  coupling  and  comparison  to  in  situ 
observational  data. 

The  feedback  of  the  wave  component  to  the  ocean  component 
(i.e.,  SWAN  to  NCOM)  includes  the  SDC,  the  wave  radiation  stress 
gradients,  and  the  characteristic  velocity  and  frequency  of  the 
wave  orbital  motion  near  the  bottom.  The  wave  motion  near  the 
bottom  is  used  to  enhance  the  bottom  drag  in  shallow  water  using 
the  parameterization  described  by  Signed  et  al.  (1990)  and  Davies 
and  Lawrence  (1994).  The  wave-radiation  stress  gradients  from 
SWAN  are  applied  in  NCOM  as  a  surface  stress.  The  SDC  from 
SWAN  is  included  within  the  Coriolis  term  in  NCOM’s  momentum 
equations  (this  is  referred  to  as  the  Stokes-Coriolis  term),  is  used  to 
advect  all  the  ocean  model  fields,  and  is  included  within  NCOM’s 
continuity  equation  (these  SDC  terms  are  implemented  as  in  Ben- 
nis  et  al.,  201 1 ).  The  SDC  is  also  used  in  the  parameterization  of  the 
enhancement  of  vertical  mixing  by  LT  as  described  by  Kantha  and 
Clayson  (2004),  i.e.,  additional  shear-production  terms,  which  con¬ 
sist  of  the  product  of  the  vertical  shear  of  the  SDC  and  the  vertical 
turbulent  momentum  flux  (which  are  referred  to  as  Stokes  produc¬ 
tion  terms),  are  added  to  the  turbulent  kinetic  energy  (TKE)  and 
vertical  turbulent  length-scale  (TLS)  equations  in  the  Mellor-Yam- 
ada  Level  2.5  (MYL2.5)  turbulence  model  used  to  compute  the  ver¬ 
tical  mixing  coefficients  in  NCOM  (see  Appendix  B). 

The  primary  focus  of  the  wave-to-ocean  coupling  in  this  study 
is  on  the  effects  of  the  SDC  in  TC  conditions.  Although  the  physical 
effects  of  the  SDC  had  not  been  thoroughly  studied  until  recently, 
early  studies  assumed  the  SDC  to  be  primarily  aligned  with  the 
surface  wind  stress.  However,  Hanley  et  al.  (2010)  found  that  seas 
are  rarely  in  wind-wave  equilibrium  and  always  produce  some 
misalignment  between  the  wind  and  the  waves.  In  particular, 
propagating  swell  waves  are  susceptible  to  misalignment  with 
the  wind  direction  (Hanley  et  al.,  2011).  Van  Roekel  et  al.  (2012) 
show  in  their  idealized  large-eddy  simulation  (LES)  studies  that  a 
greater  misalignment  angle  between  the  wind  and  the  waves  de¬ 
creases  the  generation  of  LT  and  that,  in  extreme  cases,  when  the 
SDC  directly  opposes  the  wind  stress,  there  may  be  minimal  gener¬ 
ation  of  LT  and  an  associated  decrease  in  upper-ocean  mixing.  In 
terms  of  the  sheer  structure  of  a  TC,  large  wind  waves  and  swell 
combined  with  quick  changes  in  the  wind  direction  associated 
with  the  eye,  or  inner  core,  of  the  TC  may  produce  a  large  misalign¬ 
ment  angle  between  the  wind  stress  and  the  SDC.  The  numerical 
modeling  studies  of  Sullivan  et  al.  (2012),  which  included  LES  of 
LT,  indicate  that  misalignment  between  the  wind  and  waves  near 
the  eye  of  a  hurricane  can  result  in  reduced  mixing.  In  the  simula¬ 
tions  of  Hurricane  Ivan  conducted  here,  the  effect  of  misalignment 
of  the  wind  and  waves  on  LT  and,  consequently,  vertical  mixing  is 
parameterized  by  the  Stokes  TKE  production  terms  in  the  MYL2.5 
turbulence  model  (Kantha  and  Clayson,  2004)  as  discussed  in 
Appendix  B.  Calculations  of  the  COAMPS-TC  SDC  and  subsequent 
TKE  generation  associated  with  the  SDC  are  presented  in  Section  4. 


3.  COAMPS-TC  model  configuration  and  observational  data 

3.1.  Model  setup 

The  coupled  COAMPS-TC  system  consists  of  two  independent, 
three-dimensional,  variational  (3DVAR)  atmosphere  and  ocean 
data  assimilation  systems,  three  forecast  models  (atmosphere, 
ocean,  and  wave),  and  a  coupler  that  functions  as  a  router  to  dis¬ 
tribute  the  model  forecast  fields  between  each  pair  of  forecast 
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components.  Details  of  the  model  operation  and  the  associated  air- 
ocean  exchanges  are  discussed  in  Chen  et  al.  (2010).  The  atmo¬ 
spheric  data  assimilation  is  a  special  version  of  the  Navy  Atmo¬ 
spheric  Variational  Data  Assimilation  System  (NAVDAS,  Daley 
and  Barker,  2001)  that  includes  TC  relocation  and  initialization 
algorithms.  At  the  analysis  time,  the  first-guess  fields  are  relocated 
to  the  observed  TC  location.  Next,  the  new  first-guess  fields  are 
combined  with  49  synthetic  atmospheric  temperature,  moisture, 
and  wind  profiles  at  nine  vertical  levels  to  better  define  the  TC 
and  its  immediate  environment.  In  addition,  NAVDAS  uses  a  re¬ 
duced  correlation  length  scale  and  relaxation  of  the  geostrophic 
balance  near  the  TC  center  to  improve  the  covariance  used  in  the 
analysis  (Liou  and  Sashegyi,  2011).  For  NCOM,  the  Navy  Coupled 
Ocean  Data  Assimilation  (NCODA)  3DVAR  system  is  an  equivalent 
oceanographic  implementation  of  NAVDAS  (Cummings,  2005; 
Smith  et  al.,  2011).  In  two-dimensional  mode,  NCODA  analyzes 
the  SST,  sea  ice,  significant  wave  height  (SWH),  and  sea-surface 
height  (SSH)  using  a  univariate  assimilation  method.  The  three- 
dimensional  (3D)  analysis  of  sea  temperature,  salinity,  and  u 
(east/west)  and  v  (north/south)  currents  is  multivariate.  NCODA 
has  options  to  perform  assimilation  on  either  the  native  ocean 
model  vertical  grid  or  on  constant-depth  levels.  When  the 
COAMPS-TC  system  is  run  in  cycling  mode,  the  NAVDAS  and  NCO¬ 
DA  data  assimilations  are  run  sequentially  at  the  beginning  of 
every  12-h  analysis/forecast  cycle. 

The  COAMPS-TC  model  setup  for  Hurricane  Ivan  (Fig.  2)  consists 
of  a  triple-nested  atmospheric  domain  with  18-,  6-,  and  2-km  hor¬ 
izontal  resolution  and  a  total  of  60  terrain-following  vertical  levels. 
The  outer  nest  extends  from  the  equator  to  38°  N  and  from  65°  W 
to  108°  W,  with  horizontal  dimensions  of  250  x  250.  The  two  inner 
nests  translate  in  tandem  with  the  cyclone’s  vortex  center.  Atmo¬ 
spheric  boundary  conditions  are  provided  by  the  Navy  Global 
Atmospheric  Prediction  System  (NOGAPS)  model  (Rosmond  et  al., 
2002).  The  TC  vortex  initialization  in  NAVDAS  is  enabled  through 
a  TC  warning  message  that  includes  location,  wind  radii  (distance 
of  a  certain  wind  speed  threshold  from  the  center  of  circulation), 
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Fig.  2.  Atmospheric  18-km  (white),  6-km  (yellow),  and  2-km  (orange)  resolution 
grids  and  ocean/wave  4/8-km  (red)  resolution  grids  used  for  Hurricane  Ivan 
simulations.  The  two  inner  atmospheric  grids  move  in  tandem  with  the  center  of 
the  TC.  The  domain  of  the  ocean  and  wave  grids  is  large  to  allow  for  appropriate 
numerical  spin-up  of  the  cyclone  over  the  Caribbean  Sea  before  entering  the  Gulf  of 
Mexico. 


and  intensity  information  provided  by  the  National  Hurricane  Cen¬ 
ter  in  Miami,  FL  and  the  Joint  Typhoon  Warning  Center  in  Hawaii. 
Each  data-assimilation  cycle  is  initiated  using  the  prior  12-h  fore¬ 
cast  as  background,  incorporating  quality-controlled  observations 
from  radiosondes,  aircraft,  satellites,  ships,  and  surface  stations. 

The  NCOM  configuration  (Fig.  2)  consists  of  one  4-km  resolution 
nest  that  encompasses  the  COM  and  the  Caribbean  Sea.  This  ocean 
grid  extends  from  10°  N  to  31°  N  and  from  67°  W  to  98°  W,  with 
horizontal  dimensions  of  800  x  600.  A  total  of  50  vertical  levels 
are  used,  with  36  sigma  layers  between  the  surface  and  190-m 
depth  (which  are  bottom-following  in  water  shallower  than 
190  m),  and  14  fixed-depth  layers  between  190  m  and  the  maxi¬ 
mum  depth  of  5500  m.  Bathymetry  is  from  the  Navy’s  two-min- 
ute-resolution  Digital  Bathymetric  Data  Base  (DBDB2).  Initial  and 
boundary  conditions  for  NCOM  are  from  a  global  version  of  NCOM 
(Barron  et  al.,  2004)  that  was  being  run  operationally  at  the  Naval 
Oceanographic  Office  at  Stennis  Space  Center,  MS  at  this  time.  The 
NCODA  3DVAR  ingests  observational  and  global  ocean  data,  includ¬ 
ing  quality-controlled  satellite,  ship,  profiler,  and  Modular  Ocean 
Data  Assimilation  System  (MODAS)  synthetic  profile  data  (Carnes 
et  al„  1996)  for  each  update  cycle.  The  ocean  model  includes  tides 
and  river  inflows. 

The  SWAN  model  configuration  consists  of  one  8-km  resolu¬ 
tion  nest,  with  horizontal  dimensions  of  400  x  300,  that  encom¬ 
passes  the  same  region  as  NCOM.  The  bathymetry  used  for 
SWAN  is  also  from  DBDB2.  SWAN  is  run  with  36  discrete  direc¬ 
tions  and  25  frequency  bands  to  resolve  the  wave  spectra.  The 
boundary  conditions  for  SWAN  would  normally  consist  of  energy 
spectra  from  a  larger,  global  model,  such  as  WAVEWATCH  HI 
(WW3).  However,  since  the  TC  generates  the  dominant  wave  con¬ 
tributions  to  the  wave  field  and  the  GOM  is  a  semi-enclosed  sea, 
boundary  conditions  for  SWAN  were  not  prescribed  in  these 
simulations. 

The  TC  intensity  is  sensitive  to  ocean  and  wave  model  feedback 
to  the  atmosphere  (Chen  et  al.,  2010;  Halliwell  et  al.,  201 1 ).  In  this 
study,  the  interactions  between  NCOM  and  SWAN  were  explored 
through  both  coupled  and  uncoupled  ocean-wave  simulations, 
while  the  ocean-  and  wave-model  feedback  to  the  atmosphere, 
as  discussed  in  Section  2,  was  retained  in  both  cases. 

Hurricane  Ivan  was  a  large  and  intense  TC  in  the  central  and 
eastern  GOM  from  0000  UTC  14  September  until  landfall  at 
approximately  0700  UTC  16  September  2004.  Spin-up  of 
COAMPS-TC  with  12-h  atmosphere  and  ocean  data  assimilation  cy¬ 
cles  commenced  at  0000  UTC  1  September.  Hurricane  Ivan’s  vortex 
was  initialized  on  0000  UTC  10  September  when  the  well-devel¬ 
oped  cyclone  was  located  over  the  Caribbean  Sea.  This  provided  a 
good  initial  ocean  and  wave  state  as  Ivan  entered  the  GOM  on  14 
September.  For  both  the  coupled  and  uncoupled  model  simula¬ 
tions,  a  72-h  forecast  of  Hurricane  Ivan  was  generated  at  0000 
UTC  and  1200  UTC  14  September  to  provide  comparison  with  the 
in  situ  ocean  and  wave  observational  data. 

3.2.  Observational  data  description 

Observational  data  from  the  GOM  were  supplied  by  a  total  of 
fourteen  ADCPs  (six  of  which  were  utilized  in  this  study),  five 
NDBC  buoys,  one  SRA  flight  track,  and  eight  satellite  altimeter 
passes.  The  combination  of  these  datasets  provided  a  comprehen¬ 
sive  spatial  and  temporal  sampling  of  Ivan  as  it  traversed  the  GOM. 

The  ADCPs  were  deployed  in  May  2004  along  the  outer  conti¬ 
nental  shelf  and  slope  in  the  northeastern  GOM  (Fig.  1).  The  six 
moorings  used  in  this  study  were  deployed  along  the  shallow  outer 
shelf  in  two  rows  of  three,  each  with  a  horizontal  spacing  of  15  km. 
The  first  row  of  moorings,  denoted  as  Ml -M3,  was  at  a  depth  of 
60  m,  while  the  second  row,  M4-M6,  was  at  a  depth  of  90  m.  They 
were  deployed  in  Trawl  Resistant  Bottom  Mounts  (TRBMs),  which 
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utilize  dome-shaped  pods  known  as  "Barny”  mounts  for  their  bar¬ 
nacle-like  shape  (Perkins  et  al.,  2000).  The  moorings  were 
equipped  with  Sea-Bird  Electronics  wave/tide  gauges  and  RD 
Instruments  Workhorse  ADCPs,  which  operated  at  300  kHz.  The 
ADCP  heads  were  situated  about  0.5  m  off  the  sea  floor  and  re¬ 
corded  current  profiles  with  2-m  vertical  resolution  with  an  accu¬ 
racy  of  0.5%  ±  0.5  cm  s_1  at  15-min  intervals. 

Several  NDBC  buoys  were  in  close  proximity  to  Hurricane  Ivan 
as  it  crossed  the  central  and  northern  COM  (Fig.  3).  Buoy  measure¬ 
ments  included  atmospheric  pressure,  wind  speed  and  direction, 
air  and  water  temperature,  wave  energy  spectra,  ocean  current 
velocity,  SWH,  and  wave  period  and  direction.  Based  on  NDBC  field 
studies,  SWH  accuracies  are  within  0.2  m  or  5%,  wave  periods  are 
within  1  s,  and  wave  directions  are  within  10°. 

The  European  Space  Agency’s  GlobWave  project  consolidates 
ocean  surface  wind  and  wave  data  from  multiple  satellite  instru¬ 
ments.  This  dataset  is  quality  controlled  and  available  through  a 
single  website  for  easy  access.  Detailed  information  can  be  found 
in  Snaith  et  al.  (2010).  The  GlobWave  altimeter  data  for  the  SWAN 
model  validation  during  14-16  September  2004  were  acquired 
from  three  satellites,  the  European  Remote  Sensing  satellite  (ERS- 
2),  Envisat,  and  the  GEOSAT  Follow-On  (GFO).  These  satellite  altim¬ 
eters  provided  several  snapshots  of  the  spatial  variations  of  Ivan’s 
wind  and  wave  fields  along  its  track  over  the  COM.  Gridded  model 
wind  and  wave  data  within  30  min  of  the  satellite  passing  times 
were  spatially  interpolated  to  positions  along  the  satellite  tracks. 
Additionally,  NASA  reconnaissance  aircraft  used  a  SRA  to  collect 
high-resolution  wave  spectra  and  SSH  data  for  several  periods  dur¬ 
ing  Hurricane  Ivan’s  lifecycle.  The  wave  spectra  from  the  SRA  flight 
in  the  southern  GOM  were  used  to  quantify  the  performance  of 


SWAN  in  extreme  wind  and  wave  conditions  within  all  four  quad¬ 
rants  of  the  cyclone. 

4.  Results 

4.1.  Ocean-to-wave  coupling  in  COAMPS-TC 

Low  errors  in  the  forecast  track  and  intensity  of  Hurricane  Ivan 
allowed  acceptable  statistical  comparisons  to  be  generated  for 
evaluating  the  ocean-wave  coupling  in  COAMPS-TC.  Both  NCOM- 
to-SWAN  coupled  and  uncoupled  COAMPS-TC  simulations  (details 
in  Section  4.1.2)  produced  wind  speed  errors  that  were  generally 
5  m  s_1  or  less,  with  best-track  errors  of  less  than  20  nautical  miles 
(37  km)  for  most  of  the  forecast  track.  However,  the  track  errors  in¬ 
creased  just  before  landfall  along  the  northern  GOM  coast  (Fig.  3). 
The  forecast  track  was  nearly  identical  for  all  the  simulations  of 
Hurricane  Ivan  as  the  cyclone  traversed  the  GOM  for  both  the 
0000  UTC  and  1200  UTC  14  September  2004  forecasts,  which  al¬ 
lowed  for  direct  comparisons  between  the  coupled  and  uncoupled 
model  simulations. 

4.1.3.  Near-surface  current  response  for  Hurricane  Ivan  in  COAMPS-TC 

Validation  of  the  near-surface  currents  in  NCOM  in  such  ex¬ 
treme  wind  conditions  is  quite  rare,  but  necessary  in  order  to  eval¬ 
uate  the  ocean-to-wave  coupling  performance  of  NCOM  and 
SWAN.  NCOM  was  allowed  to  spin  up  with  NCODA  ocean  data 
assimilation  for  several  weeks  prior  to  the  arrival  of  Hurricane  Ivan 
in  the  GOM  to  allow  pertinent  ocean  circulation  and  surface  fea¬ 
tures  to  develop,  thus  providing  a  good  initial  state  for  the  cy¬ 
clone’s  passage  across  the  GOM.  As  stated  previously,  the  1200 
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Fig.  3.  The  locations  of  NDBC  buoys  42001, 42003, 42039, 42040,  and  42041  are  denoted  by  white  triangles.  Hurricane  Ivan’s  observed  best  track  is  illustrated  by  the  dashed 
gray  line  and  6-hourly  position  markers  (white  squares)  and  the  model  track  by  the  solid  black  line.  The  background  colors  indicate  the  COAMPS-TC  SWH  at  forecast  hour  41 
from  the  1200  UTC  14  September  2004  72-h  forecast.  Wind  (white)  and  wave  (black)  direction  vectors  are  plotted  (relative  magnitude  is  not  indicated).  The  red  box  shows  the 
location  of  the  Scanning  Radar  Altimeter  (SRA)  flight  in  relation  to  the  wave  buoys. 
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UTC  14  September  2004  72-h  forecast  track  of  Hurricane  Ivan 
brought  the  cyclone  directly  over  the  14  ADCP  current  profilers  de¬ 
ployed  along  the  outer  continental  shelf  and  upper  slope  of  the 
GOM  just  south  of  the  Mississippi  and  Alabama  coasts.  The  ADCP 
data  indicated  that  the  shelf  currents  followed  Ekman  dynamics 
with  overlapping  surface  and  bottom  layers  during  Ivan’s  approach 
and  transitioned  to  a  dominant  surface  boundary  layer  as  the  wind 
stress  peaked  with  Ivan’s  passage  (Teague  et  al.,  2007).  In  addition, 
Hurricane  Ivan  generated  very  strong  surface  and  subsurface  cur¬ 
rents  on  the  shelf  and  slope.  For  example,  the  ADCP  Ml  measured 
currents  in  excess  of  200  cm  s_1  during  the  forced  stage  response, 
while  currents  on  the  slope  at  depths  >50  m  commonly  exceeded 
50  cm  s-1. 

To  compare  the  observed  ocean-current  response  near  the  sur¬ 
face  in  the  COAMPS-TC  simulations,  we  followed  Kuzmic  et  al. 
(2006)  for  the  calculation  of  the  magnitude  of  the  complex  correla¬ 
tion  coefficient  (CCC,  Eq.  (1))  and  the  angular  displacement,  or 
mean  directional  error  (MDE,  Eq.  (2)),  between  the  measured  ADCP 
and  NCOM  model  currents  in  the  NCOM-to-SWAN  coupled  simula¬ 
tion  (Kundu,  1976): 

rrr  =  (U0um  +  V0  Vm)  +  i(Uo  Vm  ~  Um  Vp)  n  ^ 

V(Uo  +  Vl)  ■  V(Ui  +  O 


MDE  =  tan  1 


(t to  Vm 
(u0t/m 


tlm  Vo) 
V0Vm) 


(2) 


where  u  and  fare  the  east-west  and  north-south,  observed  (o)  and 
modeled  (m),  demeaned  velocities,  and  the  angle  brackets  represent 
a  time  average.  The  CCC  and  MDE  were  computed  for  the  six  shal¬ 
low-water  ADCPs,  which  recorded  near-surface  ocean  currents. 
Inherent  in  the  computation,  the  CCC  accounts  for  both  the  current 
speed  and  direction  in  its  calculation  of  the  correlation  coefficient. 

Table  1  compares  the  72-h  1200  UTC  14  September  forecast 
simulated  and  observed  velocities  at  the  location  of  the  top  bin 


at  each  of  the  ADCP  moorings  M1-M6  (Fig.  1 )  that  had  a  near-sur¬ 
face  current  measurement  for  the  NCOM-to-SWAN  coupled 
COAMPS-TC  simulation  (BHcpl,  see  Section  4.1.2).  The  model 
velocities  were  taken  from  the  closest  grid  point  to  the  moorings 
and  interpolated  to  the  depth  of  the  observations.  Since  the  Ml- 
M6  ADCPs  were  shallow-depth  moorings  and  recorded  current 
velocities  near  the  surface,  these  statistical  results  are  important 
for  validating  that  NCOM  was  producing  satisfactory  surface  cur¬ 
rents  for  inclusion  in  SWAN.  To  better  gauge  the  accuracy  of  the 
near-surface  ocean  response  to  Hurricane  Ivan  predicted  by 
COAMPS-TC,  the  timing  of  the  prediction  with  respect  to  the  TC 
location  was  corrected  when  performing  comparisons  between 
the  model  simulations  and  the  observations.  It  is  important  to  note 
that  these  comparisons  are  very  sensitive  to  along-track  timing  er¬ 
rors  and  cross-track  spatial  errors,  and  the  timing  adjustments 
were  done  to  reduce  the  effect  of  temporal  errors  in  the  location 
of  the  TC  in  the  simulations.  Since  the  model  track  was  quite  good, 
the  temporal  adjustment  was  calculated  by  just  shifting  the  model 
time  along  the  predicted  track  to  better  match  the  observed  TC 
location.  Only  a  small  timing  adjustment  was  needed  to  account 
for  the  COAMPS-TC  timing  error  over  the  ADCP  array. 

For  each  shallow-water  ADCP,  current  measurements  were  ta¬ 
ken  in  intervals  of  15  min.  Within  every  15-min  interval,  each 
ADCP  recorded  several  measurements.  These  measurements  were 
then  averaged  every  15  min  and  velocities  recorded  at  that  time. 
NCOM  currents  were  then  compared  to  the  ADCP  currents 
throughout  the  72-h  1200  UTC  14  September  forecast,  which  pro¬ 
vided  the  best  track  and  intensity  forecast  of  Hurricane  Ivan  over 
the  ADCP  array  in  the  coupled  NCOM-to-SWAN  simulation  (the 
uncoupled  simulation  produced  a  nearly  identical  track). 

The  time-lag-adjusted  CCC  and  MDE  for  the  shallow  ADCPs 
(M1-M6)  compare  well  with  the  model  results.  Statistics  for  the 
period  encompassing  the  forced  ocean  response  to  Hurricane  Ivan 
indicate  that  mean  CCC  values  approached  or  exceeded  0.8  and 


Table  1 

ADCP  array  ocean  current  statistical  comparisons  for  the  NCOM-to-SWAN  coupled  simulation  (BHcpl)  for  the  72-h  1200  UTC  14  September  forecast.  The  magnitude  of  the 
complex  correlation  coefficient  (CCC),  mean  directional  error  (MDE),  and  current  speed  mean  bias  (MB)  are  computed  for  the  top  bins  in  each  of  the  ADCPs  M1-M6.  The  maximum 
current  speed  at  each  ADCP  and  from  NCOM  at  the  ADCP  locations  are  also  shown. 


ADCP# 

TOP  BIN  DEPTH  (m) 

CCC  (BHcpl) 

MDE  (deg,  BHcpl) 

SPEED  MB  (ms-1) 

MAX  SPEED  COAMPS  (ms-1) 

MAX  SPEED  OBS  (ms-1) 

Ml 

6 

0.83 

1.57 

0.046 

2.04 

2.14 

M2 

4 

0.77 

8.45 

0.082 

1.75 

1.87 

M3 

6 

0.79 

12.29 

0.110 

1.61 

1.73 

M4 

10 

0.78 

10.44 

0.121 

1.77 

1.96 

M5 

11 

0.82 

10.37 

0.071 

1.82 

1.91 

M6 

9 

0.83 

15.54 

0.145 

1.69 

1.82 

Fig.  4.  (a)  Time  series  of  NCOM  Eulerian  ocean  current  speed  (ms  1)  at  ADCP  Ml.  (b)  Time  series  of  observed  current  speed  from  ADCP  Ml  (ms  1). 


T.A.  Smith  et  al./ Ocean  Modelling  69  (2013)  181-194 


187 


Fig.  5.  The  friction  velocity,  u*.  as  a  function  of  10-m  wind,  U10,  based  on  ocean 
surface  drag  coefficient  formulations,  C,  0,  from  several  investigators  and  observa¬ 
tional  data  studies.  Several  parametric  models  based  on  observations  are  shown: 
default  SWAN,  based  on  Wu  (1980)  (blue  line);  Wu  with  drag  coefficient  capped  at 
2.5e~3  (green  dashed  line);  Wu  with  drag  coefficient  capped  at  0.06' 17, 0  (red  dot- 
dashed  line);  Jarosz  et  al.  (2007)  (cyan  dashes);  Hwang  (2011)  (magenta  line);  and 
Hwang  (2011)  modified  as  described  in  the  text  (gold  dashed  line).  Observational 
data  points  shown  are:  Powell  et  al.  (2003)  (black  circles);  and  Powell  and  Ginis 
(2006)  (blue  diamonds).  Lastly,  we  include  the  friction  velocities  predicted  by  other 
wave  model  physics  (not  used  herein)  Tolman  and  Chalikov  (1996,  green  Xs);  and 
Bidlot  et  al.  (2007)  (also  described  in  Ardhuin  et  al.  (2010)  (red  crosses).  These 
numerical  model  results  are  from  WAVEWATCH  III®  (Tolman,  2009),  and  specif¬ 
ically  correspond  to  simulations  of  Hurricane  Ivan  at  the  location  of  buoy  42040. 


MDE  values  were  15°  or  less  in  the  coupled  simulation  in  the  top¬ 
most  bin  for  almost  all  of  the  shallow  ADCPs.  Additional  statistics 
shown  in  Table  1  are  the  mean  bias  (MB)  of  the  current  speed  be¬ 
tween  NCOM  and  the  ADCPs  as  well  as  the  maximum  model  and 
ADCP  speeds  for  the  72-h  forecast  period.  The  results  indicate  that 
the  current  speed  MB  was  generally  less  than  0.15  m  s-1,  while  the 
NCOM  maximum  speeds  were  slightly  less  than,  but  comparable 
to,  the  maximum  speeds  of  the  ADCP  observations  just  below  the 
surface. 

Errors  in  the  CCC  and  MDE  for  the  shallow-water  ADCP  observa¬ 
tions  can  be  attributed  to  several  factors.  The  model  track  was 
excellent  as  Ivan  approached  the  southwestern  edge  of  the  ADCP 
array  nearest  to  Mil.  However,  from  that  point  an  almost  due 
north  model  track  in  both  simulations,  which  was  slightly  west 
of  the  observed  track,  continued  past  ADCP  Mil.  Although  simple 
temporal  adjustments  to  the  simulated  TC  location  were  made,  this 
deviation  in  track  as  Ivan  exited  the  ADCP  array  likely  contributed 
to  some  of  the  MDE  statistical  errors.  Errors  in  the  initial  conditions 
and  the  representation  of  the  ocean  current  velocities  at  the  begin¬ 
ning  of  a  forecast  are  always  critical;  however,  the  use  of  3DVAR 
ocean  assimilation  leading  up  to  the  72-h  forecast  reduced  this  er¬ 
ror  with  respect  to  the  ADCPs.  A  time  series  of  the  vertical  cross 
section  of  the  Eulerian  currents  at  ADCP  Ml  also  shows  a  good 
ocean  response  beneath  the  surface  as  shown  in  Fig.  4  (not  time- 
lag  adjusted).  Overall,  the  near-surface  ocean  current  speed  and 
direction  response  for  Hurricane  Ivan  in  COAMPS-TC  was  satisfac¬ 
tory  and  indicates  that  the  NCOM  surface  (top-layer)  currents  that 
are  being  passed  to  SWAN  were  reasonable. 


4.1.2.  SWAN  sensitivity  to  wave  dissipation,  drag  coefficient,  and 
NCOM  to  SWAN  coupling 

An  evaluation  of  the  friction  velocity,  u*,  in  high-wind  condi¬ 
tions  is  important  for  the  use  of  the  new  SWAN  wind  input  and 


Fig.  6.  Time  series  of  maximum  TC  intensity  and  maximum  SWH  are  illustrated  for  several  SWAN  sensitivity  tests.  Sensitivity  tests  include  a  combination  of  both  Komen  and 
Babanin  wind  input  and  dissipation  parameterizations  and  Wu  (1980)  and  Hwang  (2011)  surface  ocean  drag  formulations.  For  this  study,  three  sensitivity  tests,  Komen-Wu 
uncoupled  (KWunc,  red),  Babanin-Hwang  uncoupled  (BHunc,  blue),  and  Babanin-Hwang  coupled  (BHcpl,  green)  were  completed.  A  portion  of  the  1200  UTC  September  2004 
72-h  forecast  is  plotted  (landfall  occurs  at  approximately  hour  40  in  the  simulation). 
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dissipation  parameterizations  presented  by  Babanin  et  al.  (2010) 
(described  in  Appendix  A)  for  TC  simulations.  The  evaluation  of 
u *  based  on  the  10-m  wind  (t710)  requires  the  specification  of  a 
10-m  ocean  surface  drag  coefficient  (C10).  There  are  many  formula¬ 
tions  of  the  ocean  surface  drag  coefficient  reported  in  the  literature 
and  Fig.  5  presents  several  investigators’  resultant  u *  formulations 
based  on  their  expressions  for  Ci0  (e.g.,  Wu,  1980;  Powell  et  al., 
2003;  Powell  and  Ginis,  2006). 

Early  C10  estimates  from  Wu  (1980)  showed  u*  monotonically 
increasing  with  increasing  UI0  based  on  data  collected  from  33 
experiments  under  neutrally-stable,  open-ocean  conditions.  How¬ 
ever,  the  value  of  C10  for  high-wind  conditions  (greater  than 
30  m  s_1)  was  hypothesized  since  their  data  did  not  extend  to  such 
high  wind  speeds.  In  fact,  SWAN  utilizes  the  C]0  formulation  from 
Wu  (1980)  in  its  current  version.  However,  since  the  Wu  (1980) 
experiments,  it  has  been  observed  in  field  measurements  that  C,0 
displays  a  saturation  trend  in  high  winds,  with  a  possible  decrease 
in  magnitude  for  wind  speeds  greater  than  40  m  s  '.  Field  studies 
in  tropical  cyclone  conditions  (Powell  et  al.,  2003;  Jarosz  et  al., 
2007;  Sanford  et  al.,  2011)  have  shown  that  u *  has  a  tendency  to 
asymptote  to  a  u *  value  near  2  m  s_1  (Fig.  5)  and  to  even  decrease 
slightly  for  winds  greater  than  50  m  s-1.  Hwang  (2011)  formulated 
an  empirical  equation  to  capture  the  saturation  and  decaying 
behavior  of  C10  in  high-wind  conditions  based  on  observations 
from  Felizardo  and  Melville  (1995),  Powell  et  al.  (2003),  and  Jarosz 
et  al.  (2007): 

C10  =  10~4(-0.0160l/j0  +  0.69717,0  +  8.058)  (3) 

The  resultant  calculation  for  u*  based  on  the  Hwang  (201 1 )  for¬ 
mulation  for  C,o  is  represented  in  Fig.  5.  In  this  study,  a  modified 
Hwang  (2011)  formulation  for  C10  was  tested  in  SWAN  as  well  as 
the  classic  Wu  (1980)  formulation.  Since  the  rapid  decrease  of  u * 
to  0  for  17,o  >  68  m  s_1  in  the  Jarosz  et  al.  (2007)  formulation  has 
not  been  verified,  a  modified  Hwang  formulation  (Fig.  5,  gold 
hatched  line)  was  obtained  by  capping  the  value  of  u *  at  the  max¬ 
imum  value  from  Hwang  (2011)  formulation  of  about  2  m  s~]  for 
17,o  >  50  m  s  '. 

The  Babanin  wave  source  terms  and  the  modified  Hwang  drag 
coefficient  formulation  in  SWAN  were  implemented  in  COAMPS- 
TC  to  test  their  effectiveness  in  a  severe  TC  wind  event  and  to 
investigate  the  sensitivity  of  SWAN  compared  to  the  first  genera¬ 
tion  Komen  et  al.  (1984)  wave  source  terms  and  Wu  (1980)  drag 
formulation.  Fig.  6  shows  a  comparison  of  the  maximum  SWH 
for  Hurricane  Ivan  for  both  the  old  and  new  SWAN  wave  source 
terms  and  drag  coefficient.  Three  simulations  were  tested:  (1)  Ko- 
men-Wu  NCOM-SWAN  uncoupled  (KWunc),  (2)  Babanin-Hwang 
NCOM-SWAN  uncoupled  (BHunc),  and  (3)  Babanin-Hwang 
NCOM-SWAN  coupled  (BHcpl).  The  maximum  intensities  of  Hurri¬ 
cane  Ivan  for  each  of  the  aforementioned  tests  were  very  compara¬ 
ble,  while  the  forecast  tracks  were  nearly  identical  to  the  forecast 
track  shown  in  Fig.  3.  Upon  immediate  inspection,  a  5-6  m  differ¬ 
ence  in  the  maximum  SWH  between  the  KWunc  and  BHunc  simu¬ 
lations  was  present.  The  implementation  of  the  Babanin  wave 
source  terms  and,  more  importantly,  the  reduced  drag  coefficient 
formulation  produced  a  SWH  result  that  was  much  more  reason¬ 
able  when  compared  to  observations  in  high  wind  conditions. 
However,  it  is  important  to  note  that  all  the  u*  formulations  shown 
in  Fig.  5  are  similar  for  winds  less  than  20  m  s_1  and  will  produce 
comparable  results  in  SWAN  at  these  lower  wind  speeds. 

The  further  addition  of  NCOM-to-SWAN  coupling  (BHcpl)  in 
COAMPS-TC  acted  primarily  to  reduce  the  maximum  SWH 
throughout  the  forecast  track  of  Hurricane  Ivan  by  an  additional 
7-2  m  (Fig.  6)  in  the  presence  of  winds  greater  than  50  m  s_1.  It 
is  notable  that  the  locations  of  maximum  SWH  in  the  BHunc  and 
BHcpl  simulations  tend  to  be  different.  Comparison  of  the  BHcpl 


simulation  with  the  KBunc  simulation  shows  a  large  overall  de¬ 
crease  in  the  maximum  SWH  of  6-8  m. 

4.1.3.  Comparisons  to  NBDC  buoys 

The  data  from  several  NBDC  buoys  were  analyzed  during  the  cy¬ 
clone’s  passage  over  the  eastern  and  central  COM  for  comparison 
with  the  Hurricane  Ivan  forecasts.  NBDC  buoy  42040,  just  offshore 
of  the  northern  COM  coast  and  within  the  ADCP  array  (Fig.  1 ),  was 
directly  impacted  by  the  inner  core  of  Hurricane  Ivan.  Surface  cur¬ 
rents  at  the  location  of  buoy  42040  are  normally  quite  weak,  gen¬ 
erally  less  than  10  cm  s-1,  but  as  the  cyclone  passed  over  the  buoy, 
the  near-surface  currents  exceeded  2  m  s  1  as  measured  by  ADCP 
Ml  (Table  1 ).  Fig.  7a  is  a  plot  of  SWH  for  both  the  BHunc  and  BHcpl 
simulations  compared  to  the  observations  at  buoy  42040.  An 
along-track  time  lag  is  evident  in  the  plot  of  about  3-6  h.  The 


Buoy  42040 


Buoy  42001 


Fig.  7.  Comparison  of  SWH  from  the  uncoupled  (blue,  BHunc)  and  coupled  (red, 
BHcpl)  1200  UTC  14  September  72-h  forecasts  with  NDBC  buoy  observations 
(black),  (a)  Buoy  42040;  (b)  Buoy  42001. 


T.A.  Smith  et  al./ Ocean  Modelling  69  (2013)  181-194 


189 


BHunc  simulation  produced  SWHs  of  17-18  m  at  this  location, 
while  the  BHcpl  simulation  produced  maximum  SWHs  of  15- 
16  m,  a  difference  of  about  2  m.  Hence,  the  addition  of  the  strong, 
hurricane-induced  surface  currents  from  NCOM  into  SWAN  low¬ 
ered  the  SWHs  at  the  location  of  buoy  42040  to  more  reasonable 
values.  Although  not  shown,  the  SWHs  in  the  KWunc  simulation 


270.5  271  271.5  272  272.5  273  273.5  274  274.5  275  275.5 

Longitude  (°) 

Fig.  8.  Flight  track  of  the  airborne  SRA  (blue  line)  starting  at  2000  UTC  14 
September  2004  and  ending  at  0400  UTC  15  September  2004.  Full  hours  are  printed 
in  boxes  and  15  min  intervals  are  denoted  with  blue  squares.  Flight  segments 
starting  with  even  (odd)  hours  are  colored  dark  (light)  blue.  The  observed  (black) 
and  simulated  (red)  TC  tracks  are  shown.  The  circles  on  the  two  tracks  represent 
hourly  locations  of  the  storm  center. 


for  buoy  42040  were  in  excess  of  18  m,  with  maximum  heights 
approaching  20  m. 

Comparisons  with  data  from  another  NDBC  buoy,  42001,  lo¬ 
cated  outside  of  Hurricane  Ivan’s  inner  core,  also  showed  improve¬ 
ment  when  NCOM  was  coupled  to  SWAN.  Surface  ocean  currents  in 
September  2004  revealed  a  very  large  warm-core  eddy  that  had 
shed  from  the  Loop  Current  and  was  centered  just  east  of  buoy 
42001.  Pre-Ivan  surface  currents  at  buoy  42001  were  on  the  order 
of  0.5  m  s_1.  Inclusion  of  the  surface  currents  from  NCOM  in  SWAN 
in  the  BHcpl  simulation  reduced  the  predicted  SWHs  by  about  1  m 
at  buoy  42001  compared  with  the  BHunc  simulation  (Fig.  7b).  For 
both  buoys  42040  and  42001 ,  the  effect  of  including  strong  surface 
currents  in  the  SWAN  calculations,  whether  hurricane-induced  or 
environmental,  acted  to  improve  the  SWHs  in  SWAN  at  these  loca¬ 
tions.  Results  from  other  buoys  in  the  BHcpl  simulation,  which 
were  far  from  the  region  of  strongest  winds  and  in  areas  of  weak 
surface  currents,  showed  only  modest  differences  from  the  BHunc 
simulation. 

4.1.4.  Scanning  radar  altimeter  (SRA)  and  satellite  altimeter 
comparisons 

Although  useful  in  point  comparisons,  the  array  of  NDBC  buoys 
in  the  GOM  cannot  adequately  capture  the  evolving  wave  field  in 
each  quadrant  of  a  tropical  cyclone.  As  a  TC  translates  in  a  certain 
direction,  it  is  unlikely  that  fixed  buoys  would  provide  enough 
information  about  the  wave  field  in  each  quadrant  relative  to  the 
TC’s  center.  In  addition,  not  all  buoys  provide  directional  wave 
spectra,  which  further  reduce  the  quantification  of  the  wave  field. 
The  SRA,  on  the  other  hand,  is  an  aircraft  instrument  with  high  spa¬ 
tial  and  temporal  resolution  along  the  flight  tracks,  and  it  has  been 
successfully  used  to  observe  the  wave  field  in  the  vicinity  of  a  TC 
(Wright  et  al.,  2001). 

For  the  purpose  of  this  study,  the  SRA  NASA  aircraft  flight  track 
that  best  fits  the  simulation  window  began  at  2008  UTC  14 


Time  (h) 


Fig.  9.  Aircraft  SRA  (black),  SWAN  (blue),  and  track-adjusted  SWAN  (red)  wave  propagation  direction  (top,  degrees)  and  SWH  (middle,  m)  are  plotted  along  with  BHcpl  track 
error  (bottom,  nautical  miles  (nm))  from  2000  UTC  14  September  2004  to  0400  UTC  15  September  2004. 
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September  and  completed  at  0347  UTC  15  September,  a  flight 
duration  of  7  h  and  39  min  (Fig.  8).  The  output  fields  from  SWAN 
were  interpolated  to  the  locations  of  the  airborne  sensor  along 
the  flight  track  at  15-min  intervals  (e.g.,  at  2000  UTC,  2015  UTC, 
2030  UTC...)  for  the  0000  UTC  14  September  forecast.  The  first 
(last)  SRA  measurements  at  2008  UTC  (0347  UTC)  were  used  for 
comparison  with  the  SWAN  data  at  2000  UTC  (0400  UTC).  The 
SRA  data  were  not  regularly  spaced  in  time,  so  the  data  chosen 
were  those  measurements  closest  to  the  regular  15-min  intervals. 
The  maximum  time  discrepancy  was  less  than  four  minutes,  except 
for  the  first  and  last  set  of  measurements.  The  two  parameters  used 
for  the  validation  of  SWAN,  the  SWH  and  mean  wave  propagation 
direction  (MWPD),  were  calculated  from  the  directional  spectra  as 
follows: 


SWH  =  4 


\J J  J  E(co,d)dcodi9 


(4) 


MWPD  =  tan-' 


f  f  sindE(oj,i7)dojdi9' 
f  f  cosi9E(co,i9)dcodfl 


(5) 


where  £(cu,0)  is  the  variance  density  spectrum  as  a  function  of  fre¬ 
quency  (to)  and  azimuthal  angle  ( 0 ).  For  the  SRA  data,  both  of  these 
parameters  were  readily  calculated.  Additionally,  if  the  wave 
spectrum  was  bimodal,  both  the  primary  and  secondary  propaga¬ 
tion  directions  were  given. 


Longitude  (°) 


A  comparison  of  the  SRA  and  SWAN  SWH  and  MWPD  for  the 
0000  UTC  14  September  BHcpl  forecast  is  shown  in  Fig.  9.  The 
SWH  mean  error  (ME)  in  the  simulation  is  1.62  m  and  the  root- 
mean-square  error  (RMSE)  is  2.55  m;  however,  if  the  comparisons 
are  adjusted  to  take  in  account  spatial  and  temporal  track  errors, 
the  SWH  ME  is  reduced  to  1.25  m  and  the  RMSE  is  reduced  to 
1.78  m.  The  largest  differences  in  SWH  occur  primarily  during 
the  peaks  of  the  simulated  SWH.  The  model  MWPD  is  fairly  sim¬ 
ilar  to  the  observed  value  (Fig.  9),  with  a  ME  of  about  -7°,  and  a 
RMSE  of  34.2°  when  adjusted  for  spatial  and  temporal  track 
errors.  The  propagation  direction  was  captured  very  well; 
however,  the  SWAN  waves  were  propagating  to  the  left  of  the 
actual  wave  field  (a  MWPD  ME  of  -24.71°  when  adjusted  for 
track  errors);  therefore,  they  were  shifted  by  about  a  quarter  of 
a  quadrant  (22.5°). 

In  addition  to  the  flight  SRA  data,  several  satellite  altimeter 
measurements  of  winds  and  SWH  were  also  compared  to  the  mod¬ 
el  results.  The  three  satellites,  Envisat,  ERS-2,  and  CFO,  provided 
numerous  passes  that  yielded  useful  observational  information. 
Fig.  10a  shows  the  ground  passes  made  by  each  of  the  satellites 
for  14-16  September  2004.  Notably,  the  position  of  each  pass 
was  such  that  no  measurements  were  made  within  the  inner  core 
of  Ivan.  Fig.  10b  shows  statistics  for  all  806  SWH  measurements  ta¬ 
ken  by  satellite  compared  with  the  1200  UTC  14  September  72-h 
BHcpl  and  BHunc  forecasts,  which  were  adjusted  for  track  error 
for  optimal  comparison.  The  BHcpl  RMSE  was  0.974  m  compared 
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Fig.  10.  (a)  Satellite  altimeter  paths  with  observed  SWH  (Hs)  for  the  1200  UTC  14  September  2004  forecast.  The  magenta  line  represents  the  observed  track  of  Hurricane  Ivan, 
(b)  Scatter  diagrams  of  model  versus  altimeter  SWH  for  all  806  observations  for  the  BHunc  (bottom  left)  and  BHcpl  (bottom  right)  simulations.  The  correlation  coefficient 
(CC),  the  linear,  least-squares  fit  (red  dashed  line),  and  the  RMSE  are  shown  on  each  plot. 
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to  1.250  m  for  BHunc,  indicating  a  reduction  in  SWH  in  the  coupled 
simulation.  These  results  are  consistent  with  the  lowering  of  the 
SWH  in  the  sensitivity  tests  when  the  ocean  currents  are  included 
in  the  wave  model  coupling. 

4.2.  Wave-to-ocean  coupling  in  COAMPS-TC:  Stokes  Drift  Current 

The  SDC  and  wind  stress  vectors  in  COAMPS-TC  were  analyzed 
to  determine  the  model’s  ability  to  resolve  SDC  interactions  at 
scales  of  a  few  km  near  the  dynamic  inner  core  of  a  TC.  The  Stokes 
production  of  TKE  during  Hurricane  Ivan  and  its  effect  on  vertical 
mixing  was  calculated  following  Kantha  and  Clayson  (2004,  see 
Appendix  B).  Although  there  are  no  in  situ  observational  data  to 
explicitly  evaluate  the  SDC  for  Hurricane  Ivan,  recent  theory  on 
the  subject  implies  that  the  SDC  is  important  in  TC  ocean-wave 
interactions  (Sullivan  et  al.,  2012). 

Fig.  1 1  shows  the  misalignment  angle  between  the  SDC  and  the 
surface  wind  stress  vectors  at  forecast  hour  43  during  the  BHcpl 
simulation.  Hurricane  Ivan  was  translating  to  the  north  at  about 
10  kts  near  the  ADCP  array  at  this  time.  The  eyewall  of  the  hurri¬ 
cane  is  circled  in  Fig.  11.  The  misalignment  angle  between  the 
SDC  and  the  wind  stress  near  the  core  of  the  hurricane  approached 
180°  near  point  A,  which  is  near  the  eyewall  on  the  left  (west)  side 
of  the  track.  The  idealized  study  of  Van  Roekel  et  al.  (2012)  indi¬ 
cated  that  a  misalignment  angle  greater  than  90°  works  to  inhibit 
mixing  directly  caused  by  Langmuir  circulation  and,  in  the  most 
extreme  case  of  a  misalignment  angle  of  180°,  there  is  minimal 
generation  of  LT.  A  comparison  of  the  Stokes  shear  production  of 
TKE  for  the  BHcpl  simulation  at  points  A  and  B  in  Fig.  11  on  the  left 
and  right  sides  of  the  hurricane  track,  respectively,  is  shown  in  Ta¬ 
ble  2.  Calculation  of  the  Stokes  production  term  at  forecast  hours 
38-50  indicated  that  at  point  A,  the  Stokes  production  term  be¬ 
comes  very  small  at  forecast  hour  42  and  even  shows  a  small  neg¬ 
ative  value  at  forecast  hour  43  at  the  time  of  the  greatest 
misalignment  angle.  Note  that  a  negative  contribution  of  the 
Stokes  production  to  the  total  production  of  TKE  implies  that  the 
inclusion  of  the  SDC  is  acting  to  decrease  the  vertical  mixing.  The 


Table  2 

Calculations  of  the  Stokes  shear  TKE  production  term  (Eq.  (8))  for  the  BHcpl 
simulation  at  points  A  and  B  in  Fig.  11. 


Forecast 

hour 

Point  A 

(xlO-4  m2  s"3) 

Point  B 

(xlO-4  m2  s"3) 

Point  A 
misalignment 
angle  (°) 

Point  B 
misalignment 
angle  (°) 

38 

7.4 

7.4 

14.5 

9.4 

39 

8.5 

8.2 

12.1 

6.0 

40 

9.0 

7.9 

10.8 

5.4 

41 

9.0 

7.6 

13.4 

4.8 

42 

1.9 

7.6 

61.2 

0.6 

43 

-0.5 

8.9 

142.0 

3.4 

44 

6.4 

7.5 

15.3 

0.1 

45 

6.3 

7.3 

4.9 

0.1 

46 

5.8 

6.1 

4.2 

1.1 

47 

5.2 

4.9 

1.1 

0.7 

48 

4.6 

4.6 

0.3 

5.0 

49 

4.2 

2.7 

0.5 

0.7 

50 

3.0 

2.3 

2.3 

2.2 

results  shown  in  Table  2  are  similar  to  the  findings  of  Sullivan 
et  al.  (2012)  in  their  LES  simulations  of  LT  during  a  hurricane, 
i.e.,  large  misalignment  angles  between  the  wind  and  waves  can 
occur  on  the  left-hand-side  of  a  hurricane  during  the  passage  of 
the  eye,  and  the  large  misalignment  angle  can  result  in  oppositely 
signed  vertical  gradients  of  the  Eulerian  current  and  the  SDC, 
which  can  significantly  reduce  the  contribution  of  LT  to  vertical 
mixing.  Hence,  the  Kantha  and  Clayson  (2004)  parameterization 
of  the  enhancement  of  vertical  mixing  by  LT  using  Stokes  produc¬ 
tion  terms  in  the  TKE  and  TLS  equations  provides  a  behavior  of  LT 
with  respect  to  the  alignment  between  the  wind  and  the  waves 
that  is  consistent  with  recent  LES  studies. 

The  SDC  can  also  be  an  important  component  of  the  total  cur¬ 
rent  velocity  in  areas  with  significant  waves.  Fig.  12  shows  a 
time-depth  plot  of  the  SDC  in  the  upper  15  m  for  forecast  hours 
25-50  at  point  A  in  Fig.  11,  which  is  located  near  the  west  side 
of  the  eyewall  of  Hurricane  Ivan  at  forecast  hour  40.  The  near¬ 
surface  SDC  became  quite  large  as  the  eyewall  passed  directly  over 


Fig.  11.  Misalignment  angle  (degrees)  between  the  SDC  and  surface  wind  stress  vectors  at  forecast  hour  43  in  the  BHcpl  simulation.  The  Stokes  shear  TKE  production  is 
analyzed  at  points  A  and  B.  Hurricane  Ivan's  eyewall  is  denoted  by  the  circle,  and  the  translation  direction  and  speed  are  also  indicated  on  the  plot. 
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Vertical  Profile  of  Stokes'  Drift  Current  at  Point  A 


25  30  35  40  45 


Time  (hr) 

Fig.  12.  Time-depth  plot  of  SDC  velocity  computed  from  the  SWAN  spectra  at  point 
A  in  Fig.  11. 

point  A,  i.e.,  near  0.5  m  s-1,  at  which  time  the  surface  winds  ex¬ 
ceeded  50ms~\  The  maximum  near-surface  Eulerian  currents, 
as  discussed  in  Section  3.1.1  and  shown  in  Table  1,  are  approxi¬ 
mately  2ms~\  indicating  that  the  SDC  may  be  as  large  as  10- 
25%  of  the  Eulerian  current  near  the  surface  in  extreme  TC 
conditions. 

5.  Summary 

The  COAMPS-TC  model  simulations  for  Hurricane  Ivan  show¬ 
cased  the  state-of-the-art  ESMF  coupling  of  three  independent 
models:  the  COAMPS  atmospheric  model,  NCOM,  and  SWAN.  The 
unique  and  comprehensive  observational  data  set  for  Hurricane 
Ivan  allowed  for  the  evaluation  of  model  performance  based  on  re¬ 
cent  improvements  to  the  atmospheric,  oceanic,  and  wave  physics, 
while  gaining  a  general  but  improved  understanding  of  the  pri¬ 
mary  effects  of  ocean-wave  model  coupling  in  high-wind  condi¬ 
tions.  The  new  wind  input  and  dissipation  source  terms  (Babanin 
et  al.,  2010;  Rogers  et  al.,  2012)  and  wave  drag  coefficient  formula¬ 
tion  (Hwang,  2011),  based  on  field  observations,  significantly  im¬ 
proved  SWAN’s  wave  forecasts  for  the  simulations  of  Hurricane 
Ivan  conducted  in  this  study.  In  addition,  the  passing  of  ocean  cur¬ 
rent  information  from  NCOM  to  SWAN  further  improved  the  TC 
wave  field. 

Future  work  with  COAMPS-TC  will  continue  to  focus  on  the  air- 
ocean-wave  parameterizations  to  improve  TC  track  and  intensity. 
Recent  advances,  shown  here  in  COAMPS-TC,  include  the  addition 
of  a  wave  component,  SWAN  (Booij  et  al.,  1999),  to  exchange  wave- 
related  properties  such  as  the  Stokes’  Drift  Current  (SDC),  wave- 
radiation  stress  gradients,  and  the  near-bottom  wave  motions  that 
affect  bottom  drag  in  shallow  water.  Although  a  generalized  and 
simplified  SDC  and  TKE  study  is  shown  here  for  Hurricane  Ivan, 
additional  work  is  needed  to  better  understand  the  effects  of 
wave-to-ocean  model  coupling  on  the  upper  ocean  in  high  wind 
TC  conditions.  Additionally,  an  examination  of  wave-to-atmo- 
sphere  model  feedback  and  different  NCOM  ocean  mixing  schemes 
is  currently  ongoing.  Recent  data  obtained  during  the  Impact  of  Ty¬ 
phoons  on  the  Ocean  in  the  Pacific  (ITOP)  study  provide  excellent 
ocean  vertical  temperature  profiles  for  several  TCs  to  evaluate 
SDC-induced  turbulence  in  NCOM.  Planned  efforts  will  also  include 
evaluation  of  the  WW3  model,  which  has  recently  been  coupled 
with  COAMPS-TC  to  provide  an  alternative  wave  model.  This  re¬ 
search  is  a  continuation  of  an  Office  of  Naval  Research  (ONR)  spon¬ 
sored  National  Oceanographic  Partnership  Program  (NOPP)  project 


with  the  Rosenstiel  School  of  Marine  and  Atmospheric  Science 
(RSMAS)  at  the  University  of  Miami  and  the  University  of  Rhode  Is¬ 
land  to  generate  a  unified,  coupled  model  interface  for  providing  a 
common  feedback  mechanism  between  different  TC  models,  such 
as  the  Hurricane  Weather  Research  and  Forecast  (HWRF)  System 
model.  The  overall  ocean-to-wave  model  coupling  results  within 
COAMPS-TC  for  Hurricane  Ivan  are  encouraging  and  provide  a  plat¬ 
form  for  continued  research  into  the  next  generation  of  air-ocean- 
wave  coupled  models. 
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Appendix  A.  SWAN  deepwater  source  functions  discussion 

The  observation-based  wind  input  and  whitecapping  source 
terms  of  Rogers  et  al.  (2012)  are  used  in  the  SWAN  simulations  pre¬ 
sented  in  this  study.  These  formulations  (referred  to  as  the  “Baba¬ 
nin”  parameterizations  in  this  paper)  are  based  on  earlier  work  by 
Donelan  et  al.  (2006),  Babanin  and  Young  (2005),  Young  and  Bab¬ 
anin  (2006),  Tsagareli  (2009),  Tsagareli  et  al.  (2010),  and  Babanin 
et  al.  (2010).  The  wind  input  source  terms  in  the  new  parameteri¬ 
zation  were  taken  directly  from  the  observational  work  of  Donelan 
et  al.  (2006)  in  high  wind  conditions  over  Lake  George,  Australia, 
and  modified  to  scale  with  the  friction  velocity  u*.  The  new  dissipa¬ 
tion  function  is  observation-consistent  insofar  as  it  conforms  to 
two  features  of  dissipation  in  the  real  ocean  as  reported  in  the  lit¬ 
erature  during  the  past  decade.  The  first  feature  has  a  two-phase 
dissipation,  with  waves  of  any  particular  frequency  dissipating 
either  due  to  (1 )  the  instability  (and  breaking)  of  waves  of  that  fre¬ 
quency  or  (2)  the  destabilization  by  larger  breaking  waves  (e.g., 
through  turbulence).  The  second  feature  establishes  a  wave-break¬ 
ing  threshold,  such  that  when  the  local  spectral  density  falls  below 
a  spectral  threshold,  no  breaking  occurs  at  that  frequency.  These 
two  features  were  not  included  together  in  any  numerical  model 
until  very  recently  (by  Tsagareli,  2009  in  an  academic  model  and 
by  Ardhuin  et  al.,  2010  in  WW3).  Both  features  contrast  sharply 
with  the  dissipation  terms  of  the  previous  generation  (e.g.,  Komen 
et  al.,  1984;  WAMD1  Group,  1988;  Booij  et  al.,  1999),  in  which  all 
waves  were  considered  breaking  at  all  times  and  every  wave  sys¬ 
tem  affected  the  strength  of  dissipation  of  all  other  systems  in  a 
physically  implausible  manner  (Rogers  et  al.,  2003).  The  third  ma¬ 
jor  source  function  in  deep  water  is  that  for  four-wave  nonlinear 
interaction.  For  this  term,  we  use  the  Discrete  Interaction  Approx¬ 
imation  (Hasselmann  et  al.,  1985),  as  is  customary  for  operational 
wave  models  of  this  type  (WAMDI 1988;  Tolman,  2009;  Booij  et  al., 
2009). 

Appendix  B.  Parameterization  of  vertical  mixing  by  Langmuir 
circulation/turbulence 

Kantha  and  Clayson  (2004,  KC04)  parameterize  the  enhance¬ 
ment  of  vertical  mixing  by  Langmuir  circulation  and  Langmuir 
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turbulence  in  the  MYL2.5  turbulence  model  by  adding  Stokes 
shear-production  terms  to  the  TKE  and  TLS  equations.  These  terms 
consist  of  the  product  of  the  vertical  turbulent  momentum  flux  of 
the  Eulerian  current  times  the  vertical  shear  of  the  SDC  and  have 
the  form 


TKEsdc  =  -pu'W  -JL  -  pv'W  (6) 

where  u'W  and  v'W  are  the  vertical  turbulent  momentum  fluxes,  u', 
v't  and  w'  are  the  turbulent  velocity  fluctuations,  Us  and  Vs  are  the 
horizontal  components  of  the  SDC,  and  z  is  the  vertical  coordinate. 
Note  that  this  form  for  the  Stokes  production  terms  in  the  TKE  equa¬ 
tion  is  fairly  well  established  (e.g.,  McWilliams  et  al.,  1997;  Sullivan 
et  al.,  2007;  Grant  and  Belcher,  2009;  Janssen,  2012;  Sullivan  et  al., 
2012)  and  can  be  derived  (in  analogous  fashion  to  the  derivation  of 
the  standard  Eulerian  shear-production  terms  in  the  TKE  equation) 
from  Reynolds  averaging  of  the  vortex-force  and  Bernoulli-head 
terms  in  the  kinetic  energy  equations  formed  from  the  momentum 
equations  used  to  explicitly  simulate  Langmuir  circulation  (McWil¬ 
liams  et  al.,  1997). 

In  the  MYL2.5  turbulence  model,  the  vertical  turbulent  momen¬ 
tum  fluxes  are  computed  from  the  vertical  shear  of  the  Eulerian 
current 


u'W  =  -Km 


dU 

d Z 


v'W  =  -Km 


dV 

dz 


(7) 


where  Km  is  the  vertical  turbulent  mixing  coefficient  for  momentum 
computed  by  the  MYL2.5  turbulence  model  and  U  and  V  are  the  hor¬ 
izontal  components  of  the  Eulerian  velocity.  Substituting  (7)  into 
(6),  the  Stokes  production  of  TKE  is  calculated  as 


TKEsdc 


TdUdUs^  dV  9UA 
\dz  dz  +  dz  dz  ) 


(8) 


KC04  include  the  Stokes  production  terms  (8)  in  both  the  TKE 
and  TLS  equations.  In  the  TLS  equation,  KC04  scale  the  Stokes  pro¬ 
duction  terms  with  a  higher  weighting  than  the  standard  Eulerian 
shear-production  terms,  i.e„  the  Stokes  production  terms  are 
scaled  by  a  factor  of  7.2,  whereas  the  Eulerian  shear-production 
terms  are  scaled  by  a  factor  of  1.8.  Note  that  in  the  KC04  paper, 
the  value  used  to  scale  the  Stokes  production  terms  in  the  TLS 
equation  is  given  as  4,  which  is  incorrect  (Lakshmi  Kantha,  per¬ 
sonal  communication). 

The  effect  of  the  KC04  parameterization  is  that  the  presence  of 
the  SDC  tends  to  increase  the  vertical  mixing  coefficients  by 
increasing  both  the  TKE  and  the  TLS.  Note,  however,  that  if  the  ver¬ 
tical  shear  of  the  Eulerian  current  and  the  vertical  shear  of  the  SDC 
are  of  different  sign,  the  Stokes  TKE  production  terms  can  make  a 
negative  contribution  to  the  production  of  TKE.  Hence,  the  Stokes 
production  terms  are  sensitive  to  the  relative  alignment  of  the 
wind-driven,  near-surface,  Eulerian  current  and  the  SDC,  i.e„  to 
the  relative  alignment  of  the  wind  and  the  waves. 
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